New measurements of the cosmic infrared background fluctuations in 
deep Spitzer/IRAC survey data and their cosmological implications 

A. Kashlinsky 1 ' 4 ' 6 , R. G. Arendt 2 ' 4 , M. L. N. Ashby 5 , G. G. Fazio 5 , J. Mather 3 ' 4 , S. H. 
"^j". Moseley 3,4 



O 

S-H 

< 



o 

(N 



ABSTRACT 



We extend previous measurements of cosmic infrared background (CIB) fluc- 
tuations to 1 1° using new data from the Spitzer Extended Deep Survey. Two fields, 
with depths of ~ 12 hr/pixel over 3 epochs, are analyzed at 3.6 and 4.5 /mi. Maps 
of the fields were assembled using a self-calibration method uniquely suitable for 
probing faint diffuse backgrounds. Resolved sources were removed from the maps 
to a magnitude limit of mag^s ~ 25, as indicated by the level of the remaining 
shot noise. The maps were then Fourier-transformed and their power spectra were 
evaluated. Instrumental noise was estimated from the time-differenced data, and 
subtracting this isolates the spatial fluctuations of the actual sky. The power spectra 
of the source- subtracted fields remain identical (within the observational uncertain- 
ties) for the three epochs indicating that zodiacal light contributes negligibly to the 
CN 1 fluctuations. Comparing to 8 fim power spectra shows that Galactic cirrus cannot 

^ ■ account for the fluctuations. The signal appears isotropically distributed on the sky 

as required for an extragalactic origin. The CIB fluctuations continue to diverge to 
■ > 10 times those of known galaxy populations on angular scales out to < 1°. The 

low shot noise levels remaining in the diffuse maps indicate that the large scale fluc- 
tuations arise from the spatial clustering of faint sources well below the confusion 
noise. The spatial spectrum of these fluctuations is in reasonable agreement with 
an origin in populations clustered according to the standard cosmological model 
(ACDM) at epochs coinciding with the first stars era. 
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1. Introduction 

The Cosmic Infrared Background (CIB) is the collective radiation emitted throughout cos- 
mic history, including from sources inaccessible to current telescopic studies (see review by 
Kashlinsky, 2005). The latter category includes the very first luminous objects, currently a 
subject of intense investigations and great importance for astronomy and cosmology. CIB fluc- 
tuations can be more readily discerned than the actual mean level of the background allowing 
to overcome the significant Galactic and Solar system foregrounds at these wavelengths (Kash- 
linsky et al. 1996a,b, Kashlinsky & Odenwald 2000). It is generally believed now as a result of 
detailed numerical simulations of the formation of first structures in the standard cosmological 
model that the first objects to form in the Universe, the so-called Population III stars, were very 
massive stars that occupied a brief era at epochs inaccessible to direct telescopic observations 
(see review by Bromm & et al. 2009). Because distant galaxies and pre-galactic structures are 
clustered, the CIB has angular fluctuations with a distinct spectral and spatial signal including 
a possibly measurable contribution from the first stars era (Bond et al. 1986, Kashlinsky et al. 
2004, Cooray et al. 2004). The distribution on the sky of the luminous objects to form at early 
times should be considerably different from the cosmic pattern seen today, with the differences 
diverging toward large angular scales and being particularly prominent between 5' to 1°. Al- 
though the individual sources at very high z are too faint to observe on their own, fluctuations 
in the intensity of the cosmic infrared background (CIB) will reflect the distribution of those 
early objects after foreground sources are removed to sufficiently faint levels. 

In the course of the prior work we discovered significant source-subtracted CIB fluctua- 
tions on scales as large as ~ 5' in deep Spitzer IRAC (3.6-8 /im) data (Kashlinsky et al. 2005, 
2007a,b,c - hereafter KAMM1,2,3,4). Extensive details of our analysis and a multitude of tests 
it was subjected to along with a summary of requirements any reliable CIB data analysis should 
meet in the presence of significant foreground and instrumental effects is given in Arendt et al. 
(2010, hereafter - AKMM). As thoroughly documented in AKMM, we have been very careful 
with the processing of the individual IRAC exposures into complete deep mosaicked images. 
AKMM include details with the numerous thorough checks to confirm that we are not being 
misled by instrumental artifacts or variations of the zodiacal light. We test that we see simi- 
lar background fluctuations in various fields located at different ecliptic and Galactic latitudes. 
While all our background fluctuation measurements have been derived from IRAC data, the 
recent results of Matsumoto et al. (201 1) find similar background fluctuations using completely 
independent data from the Akari satellite, and with performing an independent analysis. This 
is the further confirmation that the fluctuations are not caused by any hidden error in the IRAC 
data or our analysis. 
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The residual CIB fluctuations remain after removing galaxies down to very faint levels and 
must arise from populations with a significant clustering component, but only low levels of the 
shot noise. As suggested by KAMM1 these CIB fluctuations may originate in early populations. 
Further, it was demonstrated by KAMM4 that there are no correlations between the source- 
subtracted IRAC maps and the faintest resolved sources observed with the HST ACS at optical 
wavelengths, which likely points to the high-z origin of the fluctuations, or at least to a very 
faint population not yet observed by other means. The high- z interpretation of the detected CIB 
anisotropics has received strong further confirmation in the recent Akari data analysis which, in 
addition, measured source-subtracted CIB fluctuations at 2.4 /im and pointed out that the colors 
of the fluctuations require them being produced by highly redshifted very luminous sources 
(Matsumoto et al. 2011). It is inherently important to advance the understanding of the nature 
and the epochs of the populations producing these CIB fluctuations with new studies which can 
be obtained measuring the spatial distribution of the fluctuations at much larger angular scales. 
This was one of the motivations of the SEDS program (Fazio et al. 2008) and the first results 
from it are presented below. 

This study followed the procedures outlined in AKMM and our previous publications 
(KAMM1-4) and the paper is organized as follows: Sec. [2] describes the data and its assem- 
bly, followed by Sec. [3] which lays out the quantities computed from the assembled data. In 
Sec. |4]we present the results on the spatial spectrum of the source-subtracted fluctuation in the 
diffuse light remaining after subtracting the instrument noise. Sec. [5] presents comprehensive 
tests done to isolate the various non-cosmological contributions to the results and both low- 
and high-,2 contributors are discussed in Sec. [6j where it is shown that the signal cannot be 
accounted for by the observed populations of known galaxies and requires either unknown new 
populations at low-2 made of low-luminosity systems and distributed very differently from the 
rest of the present day galaxies, or is dominated by the high-z populations at epochs associated 
with "first stars" and clustered according to the established concordance ACDM model. 



2. Data assembly and processing 

The Spitzer Space Telescope is a 0.85 m diameter telescope launched into an earth-trailing 
solar orbit in 2003 (Werner et al. 2004, Gehrz et al. 2007). For nearly 6 years, as it was cooled 
by liquid He, its three scientific instruments provided imaging and spectroscopy at wavelengths 
from 3.6 to 160 fim. In the time since the He supply was exhausted, Spitzer has continued to 
provide 3.6 and 4.5 fim imaging with its Infrared Array Camera (IRAC). The IRAC camera has 
a 5' x 5' field of view, and a pixel scale of 1.2", which slightly undersampled the instrument 
beam size of ~ 2" (FWHM) (Fazio et al. 2004a). 



-4- 



The SEDS program is designed to provide deep imaging at 3.6 and 4.5 /an over a total area 
of about 1 square degree, distributed over 5 well-studied regions (Fazio et al. 2008). The area 
covered is about ten times greater than previous Spitzer coverage at comparable depth. While 
the main use of the SEDS data sets will be the investigation of the individually detectable and 
countable galaxies, the remaining backgrounds in these data are well-suited for CIB studies, by 
virtue of their angular scale, sensitivity and observing strategies. The first 2 of the SEDS fields 
to be completed have been used in this study: the UltraDeep Survey field (UDS; Program ID 
= 61041) and the Extended Groth Strip (EGS; Program ID = 61042). The locations, sizes, and 
depths of these fields are listed in Table 1 . The fields are located at moderate to high Galactic 
latitudes to minimize the number of foreground stars and the brightness of the emission from 
interstellar medium (cirrus). These fields also lie at relatively high ecliptic latitudes, which 
helps minimize the brightness and temporal change in the zodiacal light from interplanetary 
dust. The observations of each field are carried out at three different epochs, spaced 6 months 
apart. 

For CIB fluctuations analysis in the assembled maps we have selected regions of approx- 
imately uniform exposure in both Channel 1 and 2. These covered a square region of approx- 
imately 21' on the side in the UDS field and a rectangular region of ~ 8' x 1° in the EGS 
field. 

The procedure for map assembly is described in our previous papers (KAMM1-4) with 
an extensive summary including all the tests given in (AKMM). Below we briefly revisit the 
adopted formalism and its advantages over the standard methods when it comes to removing 
instrumental artifacts without introducing spurious correlations. This is critically important if 
one were to robustly measure fluctuations signals as faint as that expected from first stars era, 
8F< 0.1 nW/m 2 /sr at arcminute scales. 

The individual 100 sec frame time exposures are processed by the standard IRC calibra- 
tion pipeline, and then additionally corrected for the "column-pulldown effect" which affects 
the detector output in columns that contain very bright stars. The individual frames are then 
further processed and combined into mosaicked images (with a 1.2" pixel scale) using the self- 
calibration procedure that we have previously employed for this work. The 3.6 and 4.5 fxm 
images are processed independently. At each wavelength, the frames are also processed in 
several different groups to provide multiple images that can be used to assess random and sys- 
tematic errors. One set of three groups is obtained by separating the data by each of the three 
epochs. Comparison of results between these three sets can provide indications of systematic 
errors induced by variations in the zodiacal light over 6 month intervals. The other set of two 
groups is obtained by separating the full sequence of frames into the alternating even and odd 
frame numbers. Comparison of results from these "A" and "B" subsets provides a good di- 
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agnostic of the random instrument noise, because the A and B subsets only differ by a mean 
interval of ~ 100s. 

We applied the least-squares self-calibration procedure described by (Fixsen et al. 2000). 
The approach formalizes the calibration procedure by describing the data with parameters that 
include both the detector characteristics and the true sky intensity. The derivation of these 
parameters via a least- squares algorithm yields an optimal solution for the calibration and the 
sky intensity. In this case our chosen model is given by 

D i = S a + F p + F q (1) 

where D l represents the raw data from a single pixel of a single frame, S a is the sky intensity 
at location a, F p is the offset for detector pixel p, and F q is a variable offset for each of the 4 
readouts (alternate vertical columns of the detector) and each frame. This assumes that the sky 
intensity (S a ) and the detector offsets (F p ) are invariant during the course of the observations. 
For a data set with fixed frame times (as our IRAC data), the detector dark current is included 
in the F p term as it is indistinguishable from an offset. For data sets with multiple frame 
times, a relatively simple extension of this data model could be applied. The variable offset 
F q can absorb time-dependent behavior of the detector, but only to the extent that it can be 
characterized with a single value per frame, or in some cases, a single value per readout per 
frame. 

While technically the CIB is the sum of all extragalactic emission, here we wish to exclude 
individually detectable objects and analyze the remaining CIB which is produced by fainter 
sources. The mosaicked images were cleaned of resolved sources in two stages. In the first 
stage, we construct and subtract an iterative "Model" of all the sources in the image. This is done 
using a variant of the CLEAN algorithm, in which our iterative loop consists of (a) identifying 
the brightest pixel in the image, and (b) subtracting a scaled version of the instrument PSF 
(including the broad low-level wings) to reduce this intensity by a set fraction (chosen as 50%). 
Subtracting only a fraction of a source at each iteration allows the Model to compensate for 
sources that are intrinsically extended in size, and for inaccuracies in the PSF. Various criteria 
for selecting the final iteration are given in AKMM; all lead to very similar results. In order to 
compare the signal from different sky locations, the final iteration can be chosen to correspond 
to a given shot-noise level. The second stage involves the construction of masks which are used 
to set areas at the locations of bright sources to zero. If the Model worked perfectly, this would 
not be necessary. However, small differences between the Model's ideal PSF and the effective 
PSF of the mosaicked images lead to relatively large amplitude residuals at the locations of 
bright stars. Therefore we construct a mask which is used to zero all pixels above an effective 
4a surface brightness. All 8 neighbors of any such pixel are also set to zero. This clipping is 
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done independently of the Model- subtraction and the procedure is iterated until no new pixels 
above the fixed threshold remain. The fraction of unmasked pixels is ~ 73% which enables 
robust FFT analysis. 




Fig. 1. — Left: Histograms of fluctuations in the EGS (blue) and UDS (red) fields for Chi and 2. Right: 
Number of independent Fourier elements that went into determining the power spectrum for each field 
and the Fourier binning specified above. 

The histogram distributions of the remaining pixel fluctuations after removal of the mod- 
eled sources and masking are shown in Fig. [TJ They are to a good accuracy Gaussian with very 
small residual skewness, S = ((5F) 3 /a 3 , parameters (\S\ < 0.05). 

In order to evaluate the random noise level of the maps, alternate calibrated frames for 
each dataset were mapped into separate "A" and "B" mosaics. The difference of the A and B 
mosaics should eliminate true celestial sources and stable instrumental effects and reflect only 
the random noise of the observations. The power spectrum of the instrument noise is then 
subtracted from the power spectrum of the assembled maps to give the measured signal of the 
remaining source-subtracted CIB fluctuations in the maps. 



3. Fourier analysis of the assembled data 
3.1. Definitions 

The maps under study are clipped and masked of the resolved sources, yielding the fluctua- 
tion field, SF(x). Its Fourier transform, A(<f) = f 5F(x) exp(—ix-q)d 2 x is calculated using the 
FFT. The power spectrum is P 2 (<?) = (| A(q) | 2 ), with the average taken over all the independent 
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Fourier elements corresponding to the given q. A typical rms flux fluctuation is ^ q 2 P 2 {q) /2ti 
on the angular scale of wavelength 27r/g. The correlation function, C{9) = (5F(x)-5F(x + 9)), 
is uniquely related to P2(q) via Fourier transformation. If the fraction of masked pixels in the 
maps is too high (e.g. > 40% for IRAC maps as detailed in KAMM1-4/AKMM), the large-scale 
map properties cannot be computed using the Fourier transform and instead the maps must be 
analyzed by direct calculation of C(6), which is immune to mask effects (Kashlinsky 2007). 

Several quantities are used in computations below. The power spectrum in a single band n 
is defined as P n {q). It is computed after averaging all independent Fourier elements which lie 
inside the radial interval [q, q + dq). Since the flux is a real quantity, only one half of the Fourier 
plane is independent, so that at each q there are N q /2 independent measurements of A out of 
a full ring with N g data. Because we will compare two fields of different configurations, the 
Fourier space of each was binned at the same grid of q, the mean power was then evaluated and 

the relative (Poissonian) errors on that determination were computed as J \N q . 

After demonstrating that different fields have statistically similar power spectrum, the over- 
all power at each q was averaged over the different fields with each field weighted by the relative 
statistical uncertainty defined by its configuration (see below), i.e. P(q) = ^ PiN*/ ^ N q . 
This is equivalent to simply averaging over all | A(q) | 2 available from each field i. 

We characterize the similarity (or not) of the fluctuations measured in different chan- 
nels, or at different Epochs of observations. A quantity of further interest in this context is 
the cross-power, which is the Fourier transform of the cross -correlation function C mn (9) = 
(SF m (x) ■ 5F n (x + 9)}. The cross-power spectrum is then given by P mn (q) = (A m (g)A*(g)) = 
1Z m (q)1Z n (q) + X m (q)l n (q) with 71,1 standing for the real, imaginary parts. Note the cross- 
power of a real quantity, such as the flux fluctuation, is always real, but unlike the single (auto-) 
power spectrum the cross-power can be both positive and negative. The presence (or not) of 
the same population at two channels, (m,n), can then be probed via coherence defined as 
C nm = P "p^q)p n ™J)f • If Cmn — 1 then at the two channels the same populations produce the 
diffuse signal and vice versa. 



3.2. Power spectra of the assembled fields 

The clipped and cleaned maps were Fourier transformed and power spectra evaluated. 
Visual and numerical inspection of the maps did not show any significant presence of artifacts, 
stray light, muxbleed etc. In order to average over the different fields the Fourier space was 
binned at the same set of central wavenumbers, q, in both configurations. The values of q chosen 
for the final P(q) evaluations are shown in the right panel of Fig. [T] which displays the number 
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of independent Fourier elements which is available for power spectrum evaluations. Since the 
relative error on the determined P(q) is [|iV 9 ] _ 2, the power spectrum is not determined highly 
accurately at the outer 2 bins of the EGS field; however, when combined with the independent 
determination in the UDS field, the accuracy of the averaged P(q) is acceptable (< 30 — 40%) at 
scales < 1, 000". Only the largest scale of ~ 1° probed with the EGS field alone has ^N q = 2 
and is not measured well. 

The instrument noise was evaluated from the \{A — B) maps, where A and B correspond 
to the sequences of alternating odd/even AORs. 



1 




2-rr/q arcsec 27i/q arcsec 



Fig. 2. — Fluctuations from the (^4 + B) (black triangles) and ^(A — B) maps (red error bars). 

Figure |2] shows the total power remaining in the source-subtracted maps and the noise. 
At small scales (within the beam at ~ 3' — 5') the noise contributes significantly, as it should, 
but there is a clear excess in the fluctuation power at larger scales with the noise becoming 
progressively smaller with its contribution reducing to negligible at scales greater than ~ 10' — 
20'. 
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4. Diffuse source-subtracted fluctuations 

The resulting power spectra of the residual CIB after the above subtraction and masking 
of resolved sources is shown in Fig. [3l Here we have subtracted the instrumental noise power 
by deducting the power spectra of the differences of the A and B mosaics [|(A-B)], which 
are displayed in Fig. [21 The instrumental noise power is comparable to the the CIB power at 
scales < 3", but is much smaller than the power at larger angular scales. At both 3.6 and 4.5 
/im, the power spectra are in good agreement for EGS and UDS fields, except perhaps at the 
largest angular scale (> 1000") where the small number of effective measurements leads to 
large uncertainties on the power in each individual field. 
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Fig. 3. — The spatial spectra of mean squared fluctuations for each of the fields in this study. Red/blue 
error bars correspond to UDS/EGS fields. We verify that after subtraction of the Model, the masked 
images exhibit no correlations with the removed Model. Histograms of the unmasked pixels intensities 
are approximately Gaussian, and are shown in Fig. [TJ Power spectra P(q) of the CIB are calculated 
as the amplitudes of the Fourier transforms of the masked images. The power spectra are average in 
fixed intervals of angular frequency q to reduce uncertainties at the largest angular scales, and to enable 
direct comparison of power spectra between the two fields with different sizes and shapes. The relative 
uncertainties resulting from sample (cosmic) variance were evaluated as ^/N q /2 where N q is the number 
of Fourier elements averaged with each interval in q (Fig. [J). 

In order to check for the presence of the same source-subtracted CIB at both 3.6 and 4.5 
fim, we have computed the cross-correlation power spectrum between the two channels. The 
right panel of the figure shows the resultant cross-correlation power spectra. Their similarity to 
the full power spectra confirm that the fluctuations are correlated in wavelength for both fields 
at all angular scales. Similar amplitudes would not be expected if the signal was dominated by 
noise, or by instrumental artifacts that occur independently at each wavelength. The correlation 
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of many instrumental effects are mitigated by the fact that IRAC's 3.6 and 4.5 /mi channels use 
separate optical systems and detectors to simultaneously observe fields separated by ~ 6'. 



5. Non-cosmological contributions 

The results in Fig. [3] are in excellent agreement with our earlier measurements in five addi- 
tional fields (KAMM1-2) as shown in Appendix A, but are now measured to significantly more 
accurately and extending to much greater angular scales. The isotropy of the signal (i.e. the 
consistency of results in different fields) by itself suggests a cosmological origin of the fluctua- 
tions. Nevertheless, we describe below a multitude of additional tests that we have performed 
in order to verify that non-cosmological sources cannot explain the measured signal. 



5.1. Cross-correlating different epochs of observation 

The data were analyzed separately after combining the AORs in each of the 3 epochs of 
observations, E1-E3. Epochs El and E2 (and E2 and E3) are separated by ~ 6 months and 
the detector orientation changes by ~ 180° between each pair of epochs. If the signal arises 
from the detector, rather than the sky, the fluctuations in each of the epochs should not correlate. 
Additionally, if zodiacal light fluctuations contribute significantly to the measured signal, the 
measured fluctuations should differ appreciably between the epochs. 

The random instrument noise is negligible at scales > 10". We cross-correlated the power 
spectra at epochs separated by 6 months when the detector orientation is rotated by 1 80° . Instru- 
ment artifacts that remain in a fixed geometry with respect to the detector, such as stray light and 
ghost images, are thus projected on to the field in very different patterns at the different epochs. 
The cross-correlations power spectra are shown in Fig. |4]and demonstrate that the signal comes 
from the sky rather than the detector. 
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Fig. 4. — Auto (m = n) and cross-correlation (m ^ n) power spectra between three epochs 
for EGS and UDS field at the two IRAC channels. Black, blue and red symbols correspond to 
El x E2, El x E3, E2 x E3 respectively. In the auto-correlation correlation panels the plus signs, 
asterisks and diamonds correspond to El, E2, E3 respectively. 
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5.2. Instrumental Stray Light 

One of the potential difficulties in measurements of backgrounds and their structure is 
stray light. Stray light will cause an apparent nonuniform increase in the background across the 
detector. The intensity and the pattern of the stray light depends on the distribution of sources 
outside of the field of view being observed. In the case of the IRAC 3.6 and 4.5 jum detectors, 
sources that lie < V from the top edge of the detector can scatter light back onto the detector 



(|Hora et al. 20 04). The geometry of the scattering leads to two concentrated regions of stray 
light in the upper left and right portions of the detectoiQ. Uniform diffuse illumination, such as 
the zodiacal light, leads to a fixed pattern on the detector. An individual bright source leaves 
a more concentrated and irregular patch of stray light, but at a location that is well defined 
given the known source position. An example of this is shown in the left and center panels 
of Fig. [51 Thus individual IRAC exposures have been masked^ to exclude regions that are 
potentially affected b y the stray light of sta rs identified by brightness and color thresholds from 



the 2MASS catalog (|Skrutskie et al. 20061) . However, there remains some question of whether 
or not the stray light from sources at fainter levels might yield some structure in the observed 
background. 

To address this question, we created new versions of the 3.6 /xm mosaic of the UDS field 
from two subsets for the data. The first version is made only with data from the upper half 
of detector, which can be affected by stray light if the data are insufficiently masked. The 
second version is made from only from data in the lower half of the detector which is largely 
unperturbed by stray light. The difference between these two mosaics would reveal if there 
were systematic differences between the upper and lower halves of the detector, whether due 
to stray light or any other instrumental effect. The difference image reveals no residual stray 
light near bright stars [Fig. [5](right)] and obvious large scale structures [Fig. [6]] that could arise 
from the integrated stray light of the many faint sources. Although small scale changes in the 
instrumental PSF between upper and lower halves of the detector are evident at the locations 
of bright sources. The power spectrum of the difference image is very similar to the power 
spectrum of the |(A-B) image used to assess the instrumental noise (Fig. [2]). Therefore we have 
verified that the observed large-scale fluctuations are not an artifact related to instrumental stray 
light. 



1 http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstramenthandbook/ 

2 http://irsa.ipac. caltech.edu/data/SPITZER/docs/dataanalysistools/tools/contributed/irac/straylight/ 



Fig. 5.— (left) A 10' x 10' section of the 3.6 /im mosaic of the UDS field [range = (-0.01,0.01) nW m~ 2 
sr -1 ]. (middle) A single 3.5 //m frame that shows very strong stray light from the bright star just outside 
the upper right corner of the frame. This image is displayed on a 10 x wider range [(-40,40) nW raT 2 
sr _1 ] to better show the typical structure of bright stray light, (right) A 10' x 10' section of the difference 
of mosaics that were made using (a) only data from the upper half of the detector, and (b) only data from 
the lower half of the detector. The range here is the same as in the left hand panel [(-4,4) nW m~ 2 sr -1 ]. 
No evident stray light artifacts remain. The brighter sources do exhibit small changes in the details of 
the PSF between the top and bottom halves of the array. These regions are masked in out analysis. 
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Fig. 6. — A 24' x 24' section of the difference of mosaics that were made using (a) only data from the 
upper half of the detector, and (b) only data from the lower half of the detector. The range here is the 
same as in Fig. [5] [(-4,4) nW m~ 2 sr -1 ]. If stray light from fainter sources and the background were 
present, it would appear as 5 horizontal bands matching the layout of the coverage of the field. 



-15- 



5.3. Zodiacal light 

Zodiacal light is the brightest source of foreground emissions at near-IR (Kelsall et al. 
1998), but appears to be extremely smooth with relative fluctuations below 0.2% (Abraham 
et al. 1997). More recently, Pyo et al. (2012) used mid-IR (7-24 /mi) data Akari data to 
establish upper limits of ~ 0.02% at the North Ecliptic Pole. Expected levels for zodiacal 
light fluctuations in our measurements have been discussed in the context of the Spitzer-based 
measurements (KAMM1,2 and AKMM) and also in Matsumoto et al. (2011) for the Akari- 
based results. Because the Spitzer orbits with the zodiacal cloud, the zodiacal light intensity 
and the direction of its gradient will vary with a roughly annual period in each field. The 
comparison of power spectra at, and cross-correlation power spectra between, the different 
epochs also serves as an empirical constraint on the influence of the zodiacal light. 

The fields in this study lie at high Ecliptic latitude well outside the Ecliptic plane where 
zodi is at its highest. In earlier studies it was already estimated by us to be much smaller 
than the detected fluctuations. The gradient due to zodiacal light would vary as the telescope 
moves around the Sun and hence zodiacal light would contribute different levels at different 
epochs. Then any change in the actual sky brightness or power spectra between three Epochs 
of observations separated by ~ 6 months can be attributed to changes in the zodiacal light. 

Fig. |U shows the cross-power spectra between the epochs and the signal in each of the 
epochs for both fields. The figure shows that the signal remains the same in each of the epochs 
and correlates remarkably well between them. This argues against substantial contribution from 
either the zodiacal light or detector systematics. 

To sum up this part of discussion we find no significant differences between epochs, 
demonstrating that zodiacal light is not a significant contributor to the observed spatial fluc- 
tuations. Furthermore, the relatively constant intensity of the fluctuations from 3.6 to 4.5 /xm is 
at odds with the measured zodiacal light spectrum which rises by a factor of 2.5 from 3.6 to 4.5 
fxm. Fig. |4] shows that the three separate Epochs of data produce the same fluctuation signal and 
they all correlate very well with each other. This directly confirms that the contribution of the 
time- varying component to the measured fluctuations at 3.6 and 4.5 /im is small and zodiacal 
fluctuations cannot account for the measurement. 



5.4. Galactic cirrus 

Galactic cirrus is another source of potential confusion in this measurement. However, its 
flux levels depend on the line of sight column density of the Galactic ISM and lead to variations 
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by factor of a few (~ 3—4) among all the seven fields probed in our studies (KAMM1 , KAMM2 
and SEDS). On the other hand, the measured fluctuation signal at 3.6 and 4.5 /mi is the same, 
within the statistical uncertainties, at all locations. This by itself argues against a significant 
contribution of the Galactic cirrus emission. 

Contributions of Galactic cirrus to the source-subtracted fluctuations at the near-IR wave- 
lengths have been addressed earlier in our Spitzer analysis (KAMM1-2 and AKMM) and in 
the Akari-based study (Matsumoto et al. 2011). All of the studies reached the conclusion that 
Galactic cirrus contributes well below the fluctuation levels measured there. In addition, in the 
Akari-based analysis it was demonstrated that the near-IR source-subtracted maps do not corre- 
late with the same maps observed at 100 /im, which is the most accurate probe of cirrus emis- 
sion (Arendt et al. 1998, Schlegel et al. 1998). In order to further verify that the Galactic cirrus 
emission cannot contribute appreciably to the signal measured at 3.6 and 4.5 /mi we proceed as 
outlined in KAMM1,2 and AKMM with some extra steps. We measure identical fluctuations 
in both fields in both IRAC bands despite the fact that the mean cirrus flux varies by factors of 
a few from the UDS to the EGS. The level of the fluctuations also agrees well with our earlier 
measurements at five additional locations out to angular scales of ~ 5' (KAMM1,KAMM2). 

As these high latitude deep survey fields are chosen in part based on having low ISM 
column densities, it is no surprise that we see no direct evidence of cirrus in these fields. Cir- 
rus emission at 3.6 and 4.5 /im is largely continuum emission from stochastically (transiently) 
heated small dust grains and polycyclic aromatic hydrocarbons (PAHs). An empirical measure- 
ment of cirrus emission needs to be made at longer wavelengths, such as 100 /mi where the 
emission arises from the bulk of the ISM dust (larger grains at cooler temperatures), or at 8 /mi 
where the emission is dominated by strong PAH emission bands that are closely related to the 
3.6 and 4.5 /mi emission. 

As in our prior studies, in order to evaluate cirrus contributions we have also examined 
IRAC 8 /im images generated from earlier observing programs carried out while Spitzer was still 
operating cryogenically. These data are taken with the same IRAC instrument, and have similar 
angular resolution to the 3.6 and 4.5 /im data. The 8 /mi data are somewhat shallower than the 
SEDS observations, as the IRAC 8 /mi channel was only in operation during earlier cryogenic 
observations of the EGS and UDS fields. The UDS field was covered to an average depth of 
0.67 hr under Spitzer program ID = 4002 1 , and the EGS filed was coved to a mean depth of 1 .44 
hr under Spitzer program ID = 8. The 8 /im UDS and EGS images were then masked with the 
masks from 3.6 and 4.5 /im and the power spectra computed. They are shown in Fig. [7] Their 
amplitudes at large scales are roughly proportional to the mean cirrus intensity predicted (by 
extrapolation of 100 /im observations) for these fields. However, this extrapolation is not precise 
because of intrinsic variation in the energy spectra of the cirrus. Therefore, these power spectra 
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Fig. 7. — Source subtracted CIB fluctuations at 8 /im for EGS and UDS fields. In the UDS field at 8/tm 
the maps contain a notable large-scale artifact which translates into a diagonal stripe in Fourier space. 
To probe the sensitivity of the resultant power spectrum to this feature we have computed the power 
spectrum independently from the two quadrants in the Fourier space: (q x < 0, q y > 0) (containing the 
stripe feature) and (q x > 0, q y > 0). The power spectra from these are shown with open diamonds and 
squares respectively and can be seen to be within the statistical uncertainties of the overall spectrum. 

should be viewed as setting a strict upper limit on the possible cirrus contribution at 8 /mi. The 
8 /mi fluctuations are also expected to contain contributions from the sources producing the 
measured signal at 3.6 and 4.5 /mi as well as that from other populations, so assuming that all 
of the 8 /mi fluctuations signal is produced by cirrus gives the most conservative assumption 
about its contribution. 

In order to estimate directly whether cirrus provides a significant contaminant of the mea- 
sured fluctuations, we assume that all of the measured 8 /tm diffuse emission fluctuation on 
large angular scales (where shot-noise is negligible) is generated by cirrus and then computed 
the cross-power spectrum expressing its value in terms of the coherence between channel n and 
Channel 4 at 8 /mi, C„4 = . If cirrus emission contributed all of the Channel 4 signal and 

provided a significant contribution to the measurements at 3.6 and 4.5 /mi the coherence would 
have to be C n4 ~ 1. However, the results of the computation are shown in Fig. [8] and indicate 
that, even if dominant at 8 /im, Galactic cirrus contributes negligibly to the measured signal at 
3.6 and 4.5 /im. The coherence of the EGS field is largely uncertain at the largest scale 3-4 data 
points, where the cross-power P n4 is sometimes negative, but as mentioned above, the levels of 
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cirrus fluctuations there are expected to be at a level similar to that at smaller scales^! and also 
the UDS field shows that cirrus contributes negligibly at the largest scales probed here. 

The energy spectrum of the cirrus drops by a factor of > 10 from 8 to 4.5 and 3.6 fim. 
Therefore, even if the amplitude of the large-scale power spectra at 8 /im is dominated by 
cirrus, the amplitudes at 3.6 and 4.5 /im will be at least ~ 100 times lower. At such levels, 
it remains possible that the cirrus could account for the large scale power spectra. Yet this is 
an unlikely scenario because calculation of the coherence expressed via cross-power spectra, 
Pu, -P24 (Fig. [8]) shows a relatively weak correlation between the 8 fim and shorter wavelength 
emission. This conclusion is in agreement with that of the Akari analysis (Matsumoto et al. 
2011) who demonstrate that the near-IR fluctuation do not correlate with the 100 /im data for 
the same location implying a negligible contribution of cirrus to the former. 
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Fig. 8.— Coherence, C n4 = ^^fej^ ■ for 1-4 ( blue ) and 2 " 4 ( red ) for EGS ( left ) and UDS ( ri 8 ht ) 
fields. Open symbols correspond to scales where P„4 < 0. The relative statistical uncertainties on the 
coherence resulting from cosmic variance are given by ^12/ N q which is shown with the error bars; this 
is valid at small C when all the power spectra can be assumed independent. With N q plotted in Fig. [TJ 
this uncertainty is of order 100% for the last three points in the EGS field. 



3 Since the cirrus power spectrum has power slope of n ~ — 2 (e.g. Kashlinsky & Odenwald 2000), one expects 
the level of cirrus flux fluctuations at < 1° to be similar to that at the smaller scales. 
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6. Cosmological implications 

Having established the likely cosmological origin of the measured fluctuations, we now 
turn to interpreting these results. 

There is a statistically significant signal on the remaining source- subtracted CIB fluctua- 
tions in the Spitzer data. It now extends to angular scales < 1° and, as demonstrated in Appendix 
A is consistent with our earlier measurements. KAMM1 have proposed that these CIB fluctu- 
ations originate at early epochs and KAMM3 have discussed the properties of the populations 
needed to reproduce the measured signal. On the other hand, the high-z interpretation has been 
challenged by Cooray et al. (2007) who suggested that faint ACS galaxies at intermediate z are 
responsible for the measured signal at 3.6 /im (no results for the 4.5 /mi band IRAC fluctuations 
are presented there) because they claim a drop in power at 3.6 /im when ACS -identified galaxies 
are blanked o\x% and by Thompson et al. (2007b) who proposed that the energy spectrum of the 
fluctuations is consistent with colors of stellar populations at z!8. AKMM (see Sees. 4.2, 6.3.3 
and Figs. 5, 38 there) have discussed both of these proposals (Cooray et al. 2007, Thompson et 
al. 2007b) in detail and showed that high-z interpretation remains the best description of all the 
previous data. 

Fig. [3] shows that, within the statistical uncertainties, the fields in this study have the same 
power spectrum of the source-subtracted fluctuations. Therefore, we averaged the two sets of 
results to obtain an overall power spectrum describing the CIB. The averages were weighted 
with the number of Fourier elements in each field as described above. The results of the CIB 
fluctuations after averaging over the individual fields are shown in Fig. [9j The figure indicates 
the presence of significant large-scale fluctuations remaining after removing/subtracting the 
resolved sources. 

The measured fluctuation spectrum is made of two components: small scales, <10", are 
dominated by the shot-noise produced by the discreteness of the remaining sources. The 
isotropy of the measured signal, which is further demonstrated in Appendix A for five addi- 
tional fields from our prior measurements, is consistent with it being of cosmological origin. At 
the same shot-noise level, the measured signal is in excellent agreement with our measurements 
at five other sky locations (KAMM1-2) at smaller angular scales (<300"). 



4 Two problems with the Cooray et al. analysis have been discussed as follows: 1) the power spectrum there is 
computed when over 70-80% of the field is lost to the mask which makes FFT-based analysis suspect. Kashlinsky 
(2007) shows that - when the data from Cooray et al. are analyzed via the correlation function, which is immune 
to masking - no such drop occurs. 2) Additionally, AKMM in Fig. 5 demonstrate significant inadequacies of the 
image assembly utilized in Cooray et al. 




Fig. 9. — Field-averaged CIB fluctuations at 3.6 (left), 4.5 /xm (middle) and the cross-correlation 
power spectrum. Black solid line is the contribution of the remaining known galaxies from Sullivan et 
al. (2007) who state that normal galaxies at Vega magnitudes from 22.5 - 26 can fit the observed large 
scale fluctuations. [That claim, which also appears in the comments to their arXiv:astro-ph/0609451 
posting, is contradicted by their own Figure 8 which shows the clear deficit compared to the KAMM 1 
measurements]. Because Sullivan et al. (2007) present their results only for 3.6 /xm sources, their model 
is displayed only in the left panel. Shaded areas show the residual fluctuations from Helgason et al. 
(2012), after reconstructing the near-IR CIB fluctuations of known galaxies at both 3.6 and 4.5 fim from 
a zoo of multiband galaxy luminosity function (LF) data. The shaded regions correspond to the high- 
and low faint end of the LF data. At 3.6 fim they are consistent with the black solid line although 
Helgason et al find, on average, slightly lower levels compared to Sullivan et al. (2007). The dashed 
line shows the shot-noise contribution: at 3.6 ^m the regression leads to Psn = 57.5 nJy-nW/m 2 /sr 
(or 4.8 x 10~ n nW 2 /m 4 /sr) and at 4.5 /xm the shot noise levels are Psn = 31.5 nJy-nW/m 2 /sr (or 
2.2 x 10~ n nW 2 /m 4 /sr). Since the shot noise can be expressed as Psn ~ SFcm{> mo) it is presented 
in both sets of units. Blue solid line corresponds to the high-z ACDM (toy)-model processed through 
the mask of each field and then averaged as described in the Sec. |B] It leads to the fiducial amplitude at 
5' of A 5 i = 0.07(0.05) nW/m 2 /sr at 3.6(4.5) fim. The thick solid red line shows the sum of the three 
components. 

This section discusses the constraints on the populations producing both the shot-noise 
and clustering components. It is organized as follows: we first address the limitations on the 
fluxes of the individual sources producing the large-scale clustering component which stem 
from the measured shot-noise levels and galaxy counts. Next we revisit, in light of the new 
results, the existing estimates of the levels of the clustering component of the CIB fluctuations 
from the remaining known galaxy populations and emphasize, again, that the known galaxy 
populations do not account for the large-scale clustering component, both its amplitude and 



-21 - 



spatial dependence. Finally, we discuss the limitations on the nature and epochs of the new 
populations implied by both the measured shot-noise levels and the clustering component of the 
source- subtracted CIB fluctuations. 



6.1. Galaxy counts and shot noise levels 

On scales greater than ~ 20" the fluctuation spectrum is dominated by a component due 
to clustering of the sources producing these CIB fluctuations. On scales from ~ 100" to ~ 
1000" the amplitude of this component remains roughly constant at (<5F c i us ) 2 = q 2 P(q) / (2ir) ~ 
5 x 10~ 3 nW 2 /m 4 /sr 2 at 3.6 /mi and ~ 2 x 10~ 3 nW 2 /m 4 /sr 2 at 4.5 /mi. The shot noise of 
the sources that produce these large scale fluctuations cannot exceed the observed shot noise 
component that dominates the fluctuations at scales < 10". The shot noise level is measured 
to be Psn ~ 60, 30 nJy-nW/m 2 /sr at 3.6, 4.5 /mi respectively. These levels can be compared 
to the shot noise derived from the observed galaxy counts in the fields to set lower limits on 
the magnitudes (upper limits on the flux densities) of the sources that remain unmodeled and 
unmasked in the background. 

Galaxy counts were derived independently from the reduction for the background fluctua- 
tion analysis. The Spitzer/IRAC SEDS data for the EGS and UDS fields were reduced indepen- 
dently to mosaic form using a combination of facility pipeline processing tools and IRACproc 
(Schuster et al. 2006). The resulting mosaics were constructed with 0.6" pixels using temporal 
outlier rejection to eliminate cosmic rays and instrumental artifacts. We performed source ex- 
traction using SExtractor (Betin & Arnouts 1996) in both fields and both bands - a total of four 
mosaics. This was done in dual-image mode, using the 3.6 /tm mosaic as the detection image 
and then carrying out photometry only at the positions of detected sources in the mosaics. To 
account for source confusion we inserted 10,000 randomly-placed point sources in the mosaics 
and then attempted to recover them using SExtractor with identical parameter settings as were 
used for the original source detection. In this way the source counts - which are dominated by 
galaxies below 16 Vega mag - could be corrected for completeness. For more details see Ashby 
et al. (2012), in preparation. 

The derived galaxy counts, dN{m) / dm, were integrated from a lower magnitude limit 
to infinity to derive the shot noise from the faint galaxies that may remain in the background, 
Psn(> m) = io _0 - 8m dN{m) /dm dm. Fig. [101 shows the shot noise levels from the 
measured counts as a function of the limiting magnitude, m . The shot-noise levels shown in 
Fig. [TO] are in excellent agreement with our earlier such determinations as shown in Fig. 1 of 
KAMM3. Note that the appropriate comparison there should be made to the shot noise levels 
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reached in KAMM1 using comparable integrations; our analysis of GOODS data (KAMM2) 
with 25 hr/pixel total integrations reached shot noise levels a factor of ~ 2 — 3 lower. Consistent 
with KAMM1-3 we conclude that the shot-noise level here implies that only sources at m 1 24 
remain in the cleaned maps to contribute to the source-subtracted CIB clustering component at 
scales ~ (20 — 2000)". This shot noise calculation will suffer from incompleteness in the galaxy 
counts. At faint magnitudes, mAB~21, the fields are confusion-limited (Fazio et al. 2004b) at 
the IRAC beam resolution of order ~ (3 — 4)". Correcting for incompleteness in this heavy 
confusion-limit leads to correction factors of 3> 10 and, as discussed in Fazio et al. (2004b), 
such corrections are suspect when corrections factors exceed ~ 2. Given this uncertainty, the 
limiting magnitude is likely to increase to m >25.5 — 26 after correcting for incompleteness 
as discussed in KAMM2. This is supported by Helgason et al (2012), who reconstruct galaxy 
counts from multi-wavelength luminosity function data and derive similar relations between the 
shot noise levels and the magnitude limits, m . 
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Fig. 10. — The cumulative shot noise of faint sources calculated from galaxy counts in the EGS 
(triangles) and UDS (diamonds) fields are shown as a function of the faintest magnitude to which resolved 
sources are excluded. The horizontal line indicates the measured shot noise as characterized by the power 
spectra at small spatial scales. The intersection implies that our analysis is excluding sources to tuab ~ 
24, but this is a lower limit on the magnitude because the counts are becoming significantly incomplete 
at itiab ^ 22. Reconstruction of the shot-noise using a variety of multiband galaxy luminosity function 
data from Helgason et al. (2012) is shown with shaded areas and is immune to confusion; they find the 
AB magnitude limit around ~ 25 for these shot noise levels. 

Are the sources contributing the clustered component just below the threshold of our 
source-removal, >m , or is the clustered component produced by sources significantly fainter 
with m ^> m ? The shot-noise can be expressed as Psn ~ SFcm(> m o) where S is the 
typical flux of the sources contributing to the measured fluctuations and producing the CIB 
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of net flux Fcib(> mo) (KAMM3). Thus this can be addressed by comparison between the 
observed level of the large scale fluctuations, 5F dus , and the mean CIB level implied by the 
shot noise, F CIB (> rn ) ~ Ps^/S(m ), where S = lCT a4 ( mo ~ 23 - 9 )yuJy is the typical flux of 
the population producing the large scale clustering component of the CIB. If both Psn and 
SF C \ US are dominated by the brightest sources that were not masked or modeled, m ~ 24, then 
8F c \ ns / Fcm(> mo) ~ 1 at 3.6, 4.5 /im all the way to at least ~degree scales. Thus, sources 
around m ~ 24 — 25 would need to be essentially completely clustered in structures (and voids) 
with scales of 100" — 1000". This strong large scale clustering is inconsistent with the observed 
clustering of brighter sources and the predicted clustering of faint normal galaxy populations. 
Since the shot noise level alone is consistent with normal galaxies at m > 24, we conclude 
that the large scale structure must arise from sources well below our threshold magnitude. 



6.2. Contributions from remaining galaxy populations 

To date there are three estimates of the remaining known galaxy populations to the mea- 
sured source-subtracted CIB fluctuations (KAMM1, Sullivan et al 2007, Helgason et al 2012). 
All of these estimates are consistent with each other and are discussed below. 

KAMM1 have estimated the upper limit on the CIB fluctuations due to the remaining 
known galaxies as follows. They estimated the effective limiting magnitude of the remain- 
ing galaxies and used it to evaluate net CIB flux expected from the known galaxy popula- 
tions using the counts of Fazio et al (2004). Then the upper limit on the residual CIB fluc- 
tuations were estimated by assuming that at all earlier epochs the remaining galaxies had the 
same clustering pattern as observed today on small scales with the 2-point correlation function 
£(r) = (r/5.5/i _1 Mpc) -17 . This upper limit on the clustering component of the remaining CIB 
fluctuations is shown in Fig. [9] after convolving with the IRAC beam as per Fig.l of KAMM1. 

In a different kind of analysis Sullivan et al (2007) have estimated the remaining CIB 
fluctuations at 3.6/mi for the KAMM1 parameters by reconstructing the counts in shallow, but 
wide-field, IRAC measurements as well as deep GOODS observations, encompassing the an- 
gular scales of KAMM1. They measured the clustering of resolved sources out to ~ 10' down 
to AB magnitude of ~ 24 and then used a halo model (Cooray & Sheth 2002) combined with 
conditional luminosity functions to estimate the CIB fluctuations at 3.6 jim from sources just 
below the threshold of the KAMM1 analysis. Their results are shown in Fig. [9] and are taken 
from Fig. 8 of Sullivan et al (2007). 

A different and substantially more extensive analysis was recently performed in Helgason 
et al (2012). They modelled the remaining CIB fluctuations using a massive compilation of 230 
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datasets of luminosity functions encompassing UV, optical and near-IR bands and spanning a 
wide range of redshifts. Using these data they were able to reconstruct empirically the evolution 
of the observed populations populating the redshift cone out to z ~ 5 at all wavelengths. The 
only uncertainties were due to the assumed extremes of the faint end of the luminosity function 
end, termed high-faint-end (HFE) and low-faint-end (LFE). The counts from these reconstructed 
populations were evaluated in the observer rest frames and were shown to be in excellent agree- 
ment with the measurements down to the faintest measured magnitude at all wavelengths from 
0.45 to 4.5 fim. The shot-noise reconstructed from such data is immune to confusion. This, 
coupled with the concordance cosmological model (ACDM), allowed to robustly evaluate the 
source- subtracted CIB fluctuations at 3.6 and 4.5 /mi due to the populations below the measured 
shot-noise levels for the HFE and LFE extremes. These are shown by shaded areas in Fig. [9j 

Fig. [9] shows the power spectrum expected for known galaxies from Sullivan et al. (2007) 
and Helgason et al. (2012) that are at fainter magnitudes than the shot noise constraint. This 
contribution is significantly smaller than the measured level of the fluctuations at scales greater 
than ~ 20" and also has a very different power spectrum. All of the estimates are mutually 
consistent and show that the existing galaxy populations remaining below the measured shot- 
noise levels cannot account for the clustering component of the CIB fluctuations measured by 
us. 



6.3. Unresolved very faint sources - high-z or low-z? 

The above discussion strongly suggests that the measured CIB fluctuations are produced 
by sources significantly fainter than the magnitude threshold of our source-removal procedure. 
In this case, it is likely that the measured shot-noise levels represent an upper limit on the shot- 
noise from these sources. No such faint populations have been observed in the present Universe 
and, in addition, if they exist at low-z these very faint cosmological sources would have to 
cluster very differently from the observed galaxy populations which would likely put them in 
conflict with the established ACDM cosmological model. 

On the other hand, such faint populations are expected at high- z epochs, but are beyond 
the sensitivity and/or resolution of current telescopic studies. These populations would then 
be clustered according to the standard ACDM model except their clustering properties would 
be biased because they form in haloes identified with rare density peaks of the (overall) linear 
density field at these early times (Kashlinsky et al. 2004, Cooray et al. 2004, Fernandez et 
al. 2011). To see if such populations can provide a reasonable fit to the measured large-scale 
clustering we model their power spectrum template as P(q) = P«n + A% ^ ^ ACDM j qd A ) 

y ' % PACDM(q5d A ) 
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where Pacdm is a high-z ACDM 3-dimensional template and A 5 > is to-be-fitted amplitude at 
the fiducial scale of 2rr/q 5 = 5' and cLa is the comoving angular diameter distance to z. We 
emphasize the "toy" nature of such a model: it assumes 1) that all the emission occurs over a 
narrow instant in time, and 2) the biasing is linear so that the 3-dimensional power spectrum 
of emitters is proportional to Pacdm- While the second assumption is fairly justifiable on the 
relevant linear scales, the first supposition is likely to be a significant over- simplification and in 
reality one needs to relate the projected angular power spectrum to the underlying 3-D one via 
the Limber equation after specifying the history and rate of flux production over a wide range of 
wavelengths. Thus within the framework of this "toy" model, the assumed templates are shown 
in Fig. IA-21 and are practically independent of z at the high redshifts. This is a consequence 
of the turn-over scale of the power spectra, which reflects the horizon at the matter-radiation 
equation in CDM-type models, subtending essentially the same angular scale at z 1. 

In detailed comparison of model and data, the effect of the masking on the power spectrum 
must be considered. More specifically, the multiplication of the source subtracted images by 
a mask (values of 1 or 0), corresponds to a convolution of the power spectrum by the Fourier 
transform of the mask. The effect would be redistributing some of the power from spatial scales 
where it is stronger to those where it is weaker. This has little impact on a shot noise (and 
any other white) spectrum, but can modify the appearance of a spectrum that rises sharply and 
has features of interest at large spatial scales. Thus the ACDM-type spectrum is affected more 
significantly than the spectrum of known galaxies. This is addressed in detail in Appendix B. 
Because the spectrum of the fluctuations is measured from a masked map, the model has to 
be processed correspondingly. More specifically, in Fourier space the Fourier transform of the 
model gets convolved with that of the mask and since the latter is known from our data, this 
transformation is unique. This processing is described in Appendix B, where Fig. IA-3l shows the 
fluctuation spectra after being processed through the mask. Convolution with the mask does not 
affect the shot-noise component, and affects the spectrum of known galaxies only marginally. 
However, the ACDM-type spectrum is affected more significantly by the transfer of power due 
to mask if it arises at (high) redshifts such that the subtended comoving scales reflect its steeply 
changing part. (Some of the large-scale power will then be transferred to small scales < 10", 
but this is unobservable since it is overwhelmed by the much larger shot-noise power on these 
scales.) The blue solid lines in Fig. [9] show the fits to the data over scales (100 — 1000)" 
with the effect of the masking included. The corresponding amplitudes of the model are = 
[0.07,0.05] nW/m 2 /sr at [3,6, 4.5] /mi respectively. The sum of all three components (shot-noise 
plus known galaxies from Sullivan et al. (2007) plus the high-;? ACDM population) is shown in 
Fig. |9]with the red solid line. 



7. Discussion 



The detected source- subtracted fluctuations appear of cosmological origin, but signifi- 
cantly exceed the signal from the observed known galaxy populations remaining in the data 
after the source-subtraction. The measured fluctuations are also in excellent agreement with 
our earlier measurements on the relevant scales and at the same level of shot-noise (KAMM1- 
2). Their spatial spectral distribution on sub-degree scales is also very different and is supportive 
of these fluctuations originating from sources located at early times and coincidental with the 
"first stars era". These sources are individually inaccessible to current telescopic studies and 
this measurement allows for unique characterization of their properties and spatial distribution. 
An alternative to this high-z interpretation would require the individual sources to be very faint 
as constrained by the observed shot-noise levels and the absence of correlations with the ACS 
sources (KAMM3), lie at low z and have clustering properties significantly different from those 
of observed galaxy populations. No such a new population has been observed, or proposed on 
theoretical grounds, but if true this would represent an important discovery in its own right. The 
high-2 interpretation could be ruled out using optical and shorter wavelength IR imaging such 
as would be available from cross-correlating our maps with the HST WFC3 and ACS images. 
The presence of correlations between our IRAC source- subtracted maps and those in the optical 
would argue for the new low- z populations. 

Implications of the high-z interpretation of the fluctuations have been discussed in KAMM3 
and we briefly reiterate them here in light of the new measurements. The typical biasing ex- 
pected for first haloes in the standard cosmology would lead to relative CIB fluctuations of at 
most ~ 10% on arcminute scales. Thus the measured fluctuations of less than ~ 0.1 nW/m 2 /sr 
would require the net CIB flux at 3.5 and 4.5 /mi of > (0.5 — 1) nW/m 2 /sr. Such fluxes are well 
below the claimed mean CIB excess levels from the IRTS and DIRBE measurements around 3 
/im (Dwek & Arendt 1998, Matsumoto et al. 2005) whose theoretical implications have been 
addressed in Santos et al. (2002), Salvaterra & Ferrara (2003) and Madau & Silk (2005) among 
others. These are well within the current limits on mean CIB from the 7-ray absorption mea- 
surements in low-2 blazars (Dwek et al. 2005, Aharonian et al. 2005), although potentially this 
contribution may be detectable in the high-z GRB spectra observed with Fermi LAT (Kashlin- 
sky 2005b, Kashlinsky & Band 2007, Gilmore 2011). These fluxes are also well within the 
uncertainties of CIB limits from HST deep observations of Thompson et al. (2007a). KAMM3 
argue that these populations, having surface density of ~ P SN /S 2 , must lie in the confusion 
noise of the present-day space-telescopes, and so care must be made when using filtering in as- 
sembling data from individual noisy exposures which may wash out these populations together 
with the unwanted noise (AKMM). Given the short cosmic time available at these epochs to 
produce the required CIB levels these sources must emit radiation much more efficiently with 
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much lower M/L than the present-day stellar populations (KAMM3) although there may be 
some admixture of lower-mass stars (Salvaterra et al. 2006). Even further, fluctuation measure- 
ments such as this do not directly require the sources of this emission to be exclusively stars 
and may contain additional contribution from the black-hole accretion emissions in the early 
Universe. 

This material is based upon work supported by the NASA ADP and HST-C18 grants. 
We thank Kari Helgason for many useful discussions concerning the contributions from the 
observed galaxy populations to the measured signal. 

Facilities: Spitzer (IRAC) 



Appendix 

A. Comparison with previous results 

Unlike Galactic or Solar System foregrounds, diffuse light fluctuations due to sources at 
cosmological distances should be isotropically distributed on the sky. The signal should then 
be similar at all fields provided the foreground sources are removed to approximately the same 
shot-noise levels. Fig. IA- II compares all the measurements obtained by us at comparable shot- 
noise levels (to within ~ 30%J5 For historical reasons, and because of its shallower exposure, 
the field studied in KAMM1 at 4.5 fim had populations removed only to much larger shot-noise 
level of Psn — 6 x 10 -11 nW 2 /m 4 /sr and hence is not shown in the right panel of the figure; 
however, its consistency with the KAMM2 measurements in the four GOODS fields at the same 
shot-noise level is shown in Fig. 1 of KAMM2. 



B. Correcting model templates for masking 

The power spectra are evaluated from the masked images. Although in our case the fraction 
of pixels lost to the mask is small (~ 26%), it is nonetheless important to evaluate its effects 
when comparing to predictions of a given input model. 



5 Table 1 of KAMM2 gives the shot-noise levels and Fig. 1, left of KAMM4 shows how Psn decrease with 
Model iteration. In these earlier analyses we used coarser grid of saved iterations than here, so we chose the 
shot-noise closest to the one reached in this study: Psn ~ (4.8, 2.2) x 1CP 11 nW 2 /m 4 /sr at (3.6, 4.5) /Ltm. 
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Fig. A-l. — Comparison with KAMM1/KAMM2 measurements at the shot-noise levels reached here: 
red and blue error bars represent the UDS and EGS results as shown in Fig. [3l black error bars correspond 
to the IOC/QSO 1700 field used in KAMM1, and green symbols show the four GOODS fields analyzed 
in KAMM2. Our earliest analysis in KAMM1 has not reached the shot-noise levels comparable to this 
study at 4.5 ^m and, hence the data are not shown for clarity. However, Fig. 1 of KAMM2 demonstrates 
the consistency of that measurement with the four additional GOODS field shown here at the shot-noise 
levels comparable to this study. 




Fig. A-2. — High-z shape of the concordance ACDM power spectrum projected to z = 9, 13, 17. 
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There are three potential contributions to the measured fluctuation spectrum coming from 
1) shot-noise from remaining sources, 2) clustering of known galaxies, and 3) from high-,2 
sources modelled with ACDM spectrum. The galaxy clustering contribution template we adopt 
per Sullivan et al. (2007) and Helgason et al. (2011). The template of the putative high-z 
component of sources with the ACDM is shown in Fig. IA-21 

All of these assumed components are measured from a map which is masked after clipping 
the resolved sources. Because the mask is known, this effect can and should be corrected for 
leading to unique transformation of the assumed template of the power spectrum. 

The measured fluctuations are multiplied by the mask template in real angular space. This 
is equivalent to convolution in Fourier space. Thus the shot-noise will remain flat (except it 
is now multiplied by the mask window instead of the beam), while the other components with 
P ^const will be convolved in a complicated way with the mask. To understand the effects 
of the mask and to obtain the templates after mask-processing we have conducted simulations 
generating many realizations of a diffuse field of appropriate geometric configuration with a 
given power spectrum, masking it and computing the resultant power spectrum and its distribu- 
tion. Fig. IA-31 shows the effects of the EGS and UDS masks on the simulated CIB spectra of 
the various components. 
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Table 1 . Analyzed SEDS Fields 



Region 


a, 5 


^Galj &Gal 




Size 


(£obs) 


/sky 




(deg) 


(deg) 


(deg) 




(hr) 




UDS 


34.50, -5.17 


169.98, -59.88 


30.41,-17.88 


21' x 21' 


13.6 


0.730 


EGS 


214.91,52.43 


95.95,59.81 


180.56, 60.00 


8' x 62' 


12.5 


0.725 



Note. — The fields are located at moderate to high Galactic latitudes to minimize the number 
of foregrounds stars and the brightness of the emission from interstellar medium (cirrus). These 
fields also lie at relatively high ecliptic latitudes, which helps minimize the brightness and 
temporal change in the zodiacal light from interplanetary dust. The observations of each field 
are carried out at three different epochs, spaced 6 months apart. 



